##### ANALYTICAL RESULTS #########----------------------------------------------

# Script to generate the main regression results presented in the paper. 
# The script generates the main regression Tables B1 and 1 and the main effect
# plots (Figure 9).

# Session info
# R version 4.1.0 (2021-05-18)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Big Sur 11.6.8

#-------------------------------Preparation-------------------------------------

# Load data
load("Data/vdem_panel_final.Rda")

# prepare data e.g. log and lags
vdem_panel_final <- vdem_panel_final %>% 
  mutate(oil_gas_log = log(oil_gas_value_nom + 0.01),
         v2x_polyarchy_lag = lag(v2x_polyarchy),
         e_polity2_lag = lag(e_polity2)) %>% 
  arrange(country_name, year) %>% 
  group_by(country_name) %>% 
  mutate(across(c(v2x_polyarchy, v2cademmob, v2caautmob,
                  gdppc_log_lag, growth_rate_3yrs_lag, oil_gas_log,
                  pop_log_lag, e_peaveduc, ethnic_frac,
                  intrastate_conflict, v2x_ex_party,
                  v2x_ex_military, v2x_ex_hereditary), ~ .x - lag(.x),
                .names = "{.col}_diff"),
         year_fac = factor(year)) %>% 
  ungroup()

#-------------Panel Regression (Figure 9a, 9b, Table B.1)---------------------------

# pooled model
mass.pooled <- fixest::feols(v2x_polyarchy ~ v2cademmob_osp_lag +
                               v2caautmob_osp_lag +
                               gdppc_log_lag + growth_rate_3yrs_lag +
                               oil_gas_log + pop_log_lag + 
                               e_peaveduc_lag + 
                               ethnic_frac + intrastate_conflict +
                               v2x_ex_party_lag +  v2x_ex_military_lag + 
                               v2x_ex_hereditary_lag,
                             data = vdem_panel_final,
                             panel.id = c("country_name", "year"),
                             vcov = "DK")

summary(mass.pooled)

# TW FEs 
mass.fixed <- fixest::feols(v2x_polyarchy ~ v2cademmob_osp_lag +
                              v2caautmob_osp_lag +
                              gdppc_log_lag + growth_rate_3yrs_lag +
                              oil_gas_log + pop_log_lag + 
                              e_peaveduc_lag + 
                              ethnic_frac + intrastate_conflict +
                              v2x_ex_party_lag +  v2x_ex_military_lag + 
                              v2x_ex_hereditary_lag |
                              country_name + year_fac,
                            data = vdem_panel_final,
                            panel.id = c("country_name", "year"),
                            vcov = "DK")

summary(mass.fixed)

# LDV model 
mass.LDV <- fixest::feols(v2x_polyarchy ~ v2x_polyarchy_lag +
                            v2cademmob_osp_lag + v2caautmob_osp_lag +
                            gdppc_log_lag + growth_rate_3yrs_lag +
                            oil_gas_log + pop_log_lag + 
                            e_peaveduc_lag + 
                            ethnic_frac + intrastate_conflict +
                            v2x_ex_party_lag +  v2x_ex_military_lag + 
                            v2x_ex_hereditary_lag |
                            country_name + year,
                          data = vdem_panel_final,
                          panel.id = c("country_name", "year"),
                          vcov = "DK")

summary(mass.LDV)

# FD model 
mass.FD <- plm::plm(v2x_polyarchy ~ + v2cademmob_osp_lag + v2caautmob_osp_lag +
                      gdppc_log_lag + growth_rate_3yrs_lag +
                      oil_gas_log + pop_log_lag + 
                      e_peaveduc_lag + 
                      ethnic_frac + intrastate_conflict +
                      v2x_ex_party_lag +  v2x_ex_military_lag + 
                      v2x_ex_hereditary_lag,
                    index = c("country_name", "year"),
                    model = "fd", data = vdem_panel_final)


# Export regression results to latex
texreg::htmlreg(file = "Output/table_B1.html", list(mass.pooled, mass.fixed, mass.LDV, mass.FD),
               booktabs = T,
               omit.coef = "(Intercept)", 
               custom.coef.names = c("Mobilization for democracy",
                                     "Mobilization for autocracy",
                                     "GDP per capita (log)",
                                     "GDP growth rate",
                                     "Gas and oil production (log)",
                                     "Population (log)",
                                     "Avg. years of education",
                                     "Ethnic fractionalization",
                                     "Intrastate conflict",
                                     "Ruling party dimension",
                                     "Military dimension",
                                     "Hereditary dimension",
                                     "Democracy (lag)"))


# extract all coefficients and SEs from the four models
all_model_coefs <- rbind(tidy(mass.pooled, conf.int = .95) %>% mutate(model = "Pooled"),
                         tidy(mass.fixed, conf.int = .95) %>% mutate(model = "Twoway FE"),
                         tidy(mass.LDV, conf.int = .95) %>% mutate(model = "Lagged DV"),
                         tidy(mass.FD, conf.int = .95) %>% mutate(model = "First Diffs")) %>% 
  filter(str_detect(term, 'autmob|demmob'))

# And plot them
ggplot(all_model_coefs, aes(x= term, y=estimate, group = model)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), 
                width = 0.25,size  = 1,
                position = position_dodge(width = 0.4), color="black") +
  geom_hline(yintercept = 0, color = "black",  lty = "dashed") +
  scale_x_discrete("", breaks = c("v2cademmob_osp_lag", "v2caautmob_osp_lag"),
                   labels = c("Mobilization for \n Democracy", "Mobilization for \n Autocracy")) +
  scale_y_continuous("Estimate") +
  geom_point(aes(shape = model), size = 2, position = position_dodge(width = 0.4), color="grey30") +
  scale_shape_manual("", breaks = c("Pooled", "Twoway FE",
                                    "First Diffs", "Lagged DV"),  values = c(15,16,17,18)) + 
  coord_flip() + theme_light() +
  theme(legend.position = "bottom")

# Save plot
ggsave("Output/Figure_9a.pdf",
       width = 5,
       height = 5,
       dpi = "retina")

# Calculate margina effects from the pooled model
pred1 <- plot_cap(
  mass.pooled,
  type = "response",
  condition = "v2cademmob_osp_lag", draw = F) %>% 
  select(predicted, conf.low, conf.high, mobint = condition1) %>% 
  mutate(type = "demmob")


pred2 <-plot_cap(
  mass.pooled,
  type = "response",
  condition = "v2caautmob_osp_lag", draw = F) %>%
  select(predicted, conf.low, conf.high, mobint = condition1)  %>% 
  mutate(type = "autmob")

#combine both marginal effects into one data frame
pred_all <- rbind(pred1, pred2)

# plot marginal effects
ggplot(pred_all, aes(x = mobint, y = predicted, group = type, lty = type)) +
  geom_ribbon(aes(ymin = conf.low,
                  ymax = conf.high), alpha = .75, fill = "grey70") +
  geom_line(color = "black") +
  scale_linetype_manual("", breaks = c("demmob","autmob"), values = c(1,2), 
                        labels = c("Mobilization for democracy", "Mobilization for autocracy")) +
  scale_y_continuous("Predicted EDI", limits = c(0.2,.6)) +
  xlab("Mobilization intensity") +
  theme_light() +
  theme(legend.position = "bottom")

# save plot
ggsave("Output/Figure_9b.pdf",
       width = 5,
       height = 5,
       dpi = "retina")


#-------------Survival models (Table 1, Figures B.5, B.11, B.12)----------------

# load BMR data
load("Data/bmr_auts_final.Rda")
load("Data/bmr_dems_final.Rda")

# declare Surv object
bs.surv_all_auts <- Surv(time = bmr_auts_final$t1,
                         time2 = bmr_auts_final$regime_years,
                         event = bmr_auts_final$dem_trans)

# Plot survival over time
ggsurvplot(
  fit = survfit(bs.surv_all_auts ~ 1, data = bmr_auts_final), 
  xlab = "Years", color = "black",
  ylab = "Overall survival probability",
  legend = "none", title = "Autocracies",
  surv.median.line = "hv")

# save plot
ggsave("Output/Figure_B12.pdf", width = 7, height = 5, dpi = "retina")

# Democratization

# Model 1, Cox with IV only
cfit1a <- coxph(bs.surv_all_auts ~ v2cademmob_lag + v2caautmob_lag +
                  frailty(cowcode),
                data=bmr_auts_final,
                method = "efron")

summary(cfit1a)


# Model 2, Cox with all controls and clustering
cfit2a <- coxph(bs.surv_all_auts ~ v2cademmob_lag + v2caautmob_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict + 
                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                  v2x_ex_hereditary_lag,
                cluster = cowcode,
                data=bmr_auts_final,
                method = "efron")

summary(cfit2a)

# Model 3, Cox with all controls and frailty for country
cfit3a <- coxph(bs.surv_all_auts ~ v2cademmob_lag + v2caautmob_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict + 
                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                  v2x_ex_hereditary_lag + frailty(cowcode),
                data=bmr_auts_final,
                method = "efron")

summary(cfit3a)

# Model 4, Cox with all controls and frailty for country for marginal effects
cfit4a <- coxph(bs.surv_all_auts ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict + 
                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                  v2x_ex_hereditary_lag + frailty(cowcode),
                data=bmr_auts_final,
                method = "efron")

summary(cfit4a)



# Marginal effects
me1a <- coxed(cfit4a, method="npsf", 
              newdata = mutate(bmr_auts_final,
                               v2cademmob_osp_lag=quantile(v2cademmob_osp_lag, .25, na.rm = T)),
              newdata2 = mutate(bmr_auts_final,
                                v2cademmob_osp_lag=quantile(v2cademmob_osp_lag, .75, na.rm = T)),
              id = bmr_auts_final$id)
summary(me1a, stat="mean")

me2a <- coxed(cfit4a, method="npsf",
              newdata = mutate(bmr_auts_final,
                               v2caautmob_osp_lag=quantile(v2caautmob_osp_lag, .25, na.rm = T)),
              newdata2 = mutate(bmr_auts_final,
                                v2caautmob_osp_lag=quantile(v2caautmob_osp_lag, .75, na.rm = T)),
              id = bmr_auts_final$id)
summary(me2a, stat="mean")

#  Democratic breakdown

# declare Surv object
bs.surv_all_dems <- Surv(time = bmr_dems_final$t1,
                         time2 = bmr_dems_final$regime_years,
                         event = bmr_dems_final$aut_trans)

# Plot survival probabilities
ggsurvplot(
  fit = survfit(bs.surv_all_dems ~ 1, data = bmr_dems_final), 
  xlab = "Years", color = "black",
  ylab = "Overall survival probability",
  legend = "none", title = "Democracies",
  surv.median.line = "hv")

ggsave("Output/Figure_B11.pdf", width = 7, height = 5, dpi = "retina")


# Model 1
cfit1b <- coxph(bs.surv_all_dems ~ v2cademmob_lag + v2caautmob_lag +
                  frailty(cowcode),
                data=bmr_dems_final,
                method = "efron")

summary(cfit1b)

cfit2b <- coxph(bs.surv_all_dems ~ v2cademmob_lag  +v2caautmob_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + 
                  intrastate_conflict + 
                  el_pres +  
                  v2xps_party_lag + 
                  v2juhcind_lag + cluster(id),
                data = bmr_dems_final,
                method = "efron")

summary(cfit2b)



cfit3b <- coxph(bs.surv_all_dems ~ v2cademmob_lag + v2caautmob_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict + 
                  el_pres +  v2xps_party_lag + 
                  v2juhcind_lag + frailty(cowcode),
                data=bmr_dems_final,
                method = "efron")

summary(cfit3b)

# For marginal effects
cfit4b <- coxph(bs.surv_all_dems ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict + 
                  el_pres +  v2xps_party_lag + 
                  v2juhcind_lag + frailty(cowcode),
                data=bmr_dems_final,
                method = "efron")

summary(cfit4b)

me2b <- coxed(cfit4b, method="npsf",
              newdata = mutate(bmr_dems_final,
                               v2caautmob_osp_lag=quantile(v2caautmob_osp_lag, .25, na.rm = T)),
              newdata2 = mutate(bmr_dems_final,
                                v2caautmob_osp_lag=quantile(v2caautmob_osp_lag, .75, na.rm = T)),
              id = bmr_dems_final$id)
summary(me2b, stat="mean")

# Model output to LaTeX
texreg::htmlreg(file = "Output/table_1.html",
  l = list(cfit1a, cfit3a, cfit1b, cfit3b),
  custom.model.names = c("Dem. transition", "Dem. transition",
                         "Dem. breakdown", "Dem. breakdown"),
  custom.coef.names = c("Mobilization for democracy",
                        "Mobilization for autocracy",
                        "GDP per capita (log)",
                        "GDP growth rate",
                        "Gas and oil production (log)",
                        "Population (log)",
                        "Avg. years of education",
                        "Ethnic fractionalization",
                        "Intrastate conflict",
                        "Ruling party dimension",
                        "Military dimension",
                        "Hereditary dimension",
                        "Elected president",
                        "Party institutionalization index",
                        "High court independence"))




